Автор: Винник Екатерина Петровна, 22.М04
Поскольку контуры и другие скачкообразные изменения яркости связаны с высокочастотными составляющими, повышение резкости изображения может быть достигнуто при помощи процедуры высокочастотной фильтрации в частотной области, которая подавляет низкочастотные составляющие и не затрагивает высокочастотную часть Фурье-преобразования. Фильтры высоких частот получаются вычитанием из 1 фильтров низких частот, и поэтому имеют аналогичные фильтрам низких частот свойства. Рассмотрим три фильтра высоких частот -- идеальный, Баттерворта. Гаусса: $$ \begin{equation} \begin{gathered} H(u, v) = \begin{cases} 0, \textrm{at } D(u, v) \leq D_0 \\ 1, \textrm{else}; \end{cases} \\ H(u, v) = \frac{(D(u, v) / D_0)^{2n}}{1 + (D(u, v) / D_0)^{2n}}; \\ H(u, v) = 1 - \exp (-D^2(u, v) / 2 D_0); \\ \end{gathered} \end{equation} $$ где $D(u, v)$ расстояние от центра изображения, $D_0$ константа.
import numpy as np
def get_fft(img, shift=False):
h = img.shape[0]
w = img.shape[1]
img = img.astype(complex)
x_arr = np.arange(h)
y_arr = np.arange(w)
x_mat = np.repeat(x_arr.reshape((h, 1)), w, 1)
y_mat = np.repeat(y_arr.reshape((1, w)), h, 0)
if shift:
exp_val = (x_mat + y_mat) * np.pi * 1j
img = np.multiply(img, np.exp(exp_val))
if len(img.shape) == 3:
img = np.swapaxes(img, 0, 2)
img_fft = np.fft.fft2(img)
else:
img_fft = np.fft.fft2(img)
return img_fft
import numpy as np
def get_ifft(img_fft, shift=False):
h = img_fft.shape[0]
w = img_fft.shape[1]
img_fft = img_fft.astype(complex)
u_arr = np.arange(h)
v_arr = np.arange(w)
u_mat = np.repeat(u_arr.reshape((h, 1)), w, 1)
v_mat = np.repeat(v_arr.reshape((1, w)), h, 0)
if shift:
exp_val = (u_mat + v_mat) * np.pi * 1j
img_fft = np.multiply(img_fft, np.exp(exp_val))
if len(img_fft.shape) == 3:
img_fft = np.swapaxes(img_fft, 0, 2)
img_ifft = np.fft.ifft2(img_fft)
else:
img_ifft = np.fft.ifft2(img_fft)
return img_ifft
import numpy as np
def apply_ideal_HP(img_fft, d0perc=.5):
h = img_fft.shape[0]
w = img_fft.shape[1]
d0 = d0perc * np.min([h, w]) / 2
u_arr = np.arange(h)
v_arr = np.arange(w)
u_mat = np.repeat(u_arr.reshape((h, 1)), w, 1)
v_mat = np.repeat(v_arr.reshape((1, w)), h, 0)
dist_mat = np.sqrt(np.power((u_mat - h//2), 2) + np.power((v_mat - w//2), 2))
filter_mat = np.zeros_like(dist_mat)
filter_mat[dist_mat > d0] = 1.
if len(img_fft.shape) == 3:
channels = img_fft.shape[-1]
for ch in range(channels):
img_fft[:, :, ch] = np.multiply(img_fft[:, :, ch], filter_mat)
else:
img_fft = np.multiply(img_fft, filter_mat)
return img_fft, filter_mat
import numpy as np
def apply_bw_HP(img_fft, d0perc=.5, n=1):
h = img_fft.shape[0]
w = img_fft.shape[1]
d0 = d0perc * np.min([h, w]) / 2
u_arr = np.arange(h)
v_arr = np.arange(w)
u_mat = np.repeat(u_arr.reshape((h, 1)), w, 1)
v_mat = np.repeat(v_arr.reshape((1, w)), h, 0)
dist_mat = np.sqrt(np.power((u_mat - h//2), 2) + np.power((v_mat - w//2), 2))
filter_mat = np.full_like(dist_mat, 1.) - 1 / (1 + np.power((dist_mat / d0), 2*n))
if len(img_fft.shape) == 3:
channels = img_fft.shape[-1]
for ch in range(channels):
img_fft[:, :, ch] = np.multiply(img_fft[:, :, ch], filter_mat)
else:
img_fft = np.multiply(img_fft, filter_mat)
return img_fft, filter_mat
import numpy as np
def apply_gauss_HP(img_fft, d0perc=.5):
h = img_fft.shape[0]
w = img_fft.shape[1]
d0 = d0perc * np.min([h, w]) / 2
u_arr = np.arange(h)
v_arr = np.arange(w)
u_mat = np.repeat(u_arr.reshape((h, 1)), w, 1)
v_mat = np.repeat(v_arr.reshape((1, w)), h, 0)
dist_mat = np.sqrt(np.power((u_mat - h//2), 2) + np.power((v_mat - w//2), 2))
filter_mat = np.full_like(dist_mat, 1.) - np.exp(-np.power(dist_mat, 2) / (2 * d0 * d0))
if len(img_fft.shape) == 3:
channels = img_fft.shape[-1]
for ch in range(channels):
img_fft[:, :, ch] = np.multiply(img_fft[:, :, ch], filter_mat)
else:
img_fft = np.multiply(img_fft, filter_mat)
return img_fft, filter_mat
Загрузим изображения для дальнейшей обработки.
import numpy as np
import matplotlib.pyplot as plt
from PIL import Image
img_1 = Image.open('./text1.jpg')
img_2 = Image.open('./text2.jpg')
# represent the images as tensors
np_1 = np.array(img_1)
np_2 = np.array(img_2)
print('Picture dimensions:')
print(np_1.shape)
print(np_2.shape)
Picture dimensions: (256, 256, 3) (256, 256, 3)
# simplify the images by summing the channels
np_1_sum = np.sum(np_1, axis=2) // 3
np_2_sum = np.sum(np_2, axis=2) // 3
Применим быстрое преобразование Фурье к изображениям:
np_1_fft = get_fft(np_1_sum, True)
np_2_fft = get_fft(np_2_sum, True)
# log of abs
np_1_fft_logabs = np.log(np.abs(np_1_fft)+1)
np_2_fft_logabs = np.log(np.abs(np_2_fft)+1)
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].imshow(np_1_sum, cmap='gray')
axs[0, 0].set_title('Original of the 1st image')
axs[0, 1].matshow(np_1_fft_logabs, cmap='gray')
axs[0, 1].set_title('Log of abs. of FFT of the 1st image')
axs[1, 0].imshow(np_2_sum, cmap='gray')
axs[1, 0].set_title('Original of the 2nd image')
axs[1, 1].matshow(np_2_fft_logabs, cmap='gray')
axs[1, 1].set_title('Log of abs. of FFT of the 2nd image')
plt.show()
Применим идеальный фильтр для различных значений $D_0$ (которое определяется как процент высоты изображения).
np_1_id1, _ = apply_ideal_HP(np_1_fft, .01)
np_1_id2, _ = apply_ideal_HP(np_1_fft, .05)
np_1_id3, _ = apply_ideal_HP(np_1_fft, .1)
np_1_id1i = np.abs(get_ifft(np_1_id1))
np_1_id2i = np.abs(get_ifft(np_1_id2))
np_1_id3i = np.abs(get_ifft(np_1_id3))
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].imshow(np_1_sum, cmap='gray')
axs[0, 0].set_title('Original of the 1st image')
axs[0, 1].imshow(np_1_id1i, cmap='gray')
axs[0, 1].set_title('Ideal highpass filter, D_0 = 1%H')
axs[1, 0].imshow(np_1_id2i, cmap='gray')
axs[1, 0].set_title('Ideal highpass filter, D_0 = 5%H')
axs[1, 1].imshow(np_1_id3i, cmap='gray')
axs[1, 1].set_title('Ideal highpass filter, D_0 = 10%H')
plt.show()
Наилучшие результаты достигнуты при $D_0 \sim 5\%$.
Действие идеального ФВЧ противоположно действию идеального ФНЧ в том смысле, что он обнуляет все частоты, попадающие внутрь круга радиуса $D_0$ , одновременно пропуская без ослабления все частоты, лежащие вне круга. Идеальные ФВЧ обладают такими же свойствами в отношении звона, как и идеальные ФНЧ, это вполне продемонстрировано на изображениях выше. Контуры на лучшем из результирующих изображений гораздо более ровные и меньше искажены, маленькие объекты отфильтрованы надлежащим образом.
Постоянный фон на всех изображениях, обработанных высокочастотными фильтрами равен нулю, поскольку высокочастотная фильтрация аналогична дифференцированию в пространственной области.
Применим фильтр Баттерворта при различных $D_0$ и различных значениях порядка $n$. Сначала рассмотрим различные значения $D_0$ при фиксированном значении порядка $n=2$.
## order 2
np_1_bw1_2, _ = apply_bw_HP(np_1_fft, .01, 2)
np_1_bw2_2, _ = apply_bw_HP(np_1_fft, .05, 2)
np_1_bw3_2, _ = apply_bw_HP(np_1_fft, .10, 2)
np_1_bw1i_2 = np.abs(get_ifft(np_1_bw1_2))
np_1_bw2i_2 = np.abs(get_ifft(np_1_bw2_2))
np_1_bw3i_2 = np.abs(get_ifft(np_1_bw3_2))
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].imshow(np_1_sum, cmap='gray')
axs[0, 0].set_title('Original of the 1st image')
axs[0, 1].imshow(np_1_bw1i_2, cmap='gray')
axs[0, 1].set_title('Butterworth highpass filter, D_0 = 1%H, n = 2')
axs[1, 0].imshow(np_1_bw2i_2, cmap='gray')
axs[1, 0].set_title('Butterworth highpass filter, D_0 = 5%H, n = 2')
axs[1, 1].imshow(np_1_bw3i_2, cmap='gray')
axs[1, 1].set_title('Butterworth highpass filter, D_0 = 10%H, n = 2')
plt.show()
Наилучшие результаты достигнуты при $D_0 \sim 5\%$. Как и в случае низкочастотных фильтров, высокочастотные фильтры Баттерворта ведут себя более гладко по сравнению с идеальными ФВЧ -- это вполне видно на изображениях выше. При $D_0$ контуры изображения подчеркнуты, вместе с тем отсутствует шум, вносимый идеальным ФВЧ. Переход к более высоким значениям частоты среза для ФВЧ Баттерворта совершается гораздо более плавно.
Рассмотрим результаты обработки изображения фильтром при различных значениях $D_0$ и фиксированном значении порядка $n=4$.
## order 4
np_1_bw1_4, _ = apply_bw_HP(np_1_fft, .01, 4)
np_1_bw2_4, _ = apply_bw_HP(np_1_fft, .05, 4)
np_1_bw3_4, _ = apply_bw_HP(np_1_fft, .10, 4)
np_1_bw1i_4 = np.abs(get_ifft(np_1_bw1_4))
np_1_bw2i_4 = np.abs(get_ifft(np_1_bw2_4))
np_1_bw3i_4 = np.abs(get_ifft(np_1_bw3_4))
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].imshow(np_1_sum, cmap='gray')
axs[0, 0].set_title('Original of the 1st image')
axs[0, 1].imshow(np_1_bw1i_4, cmap='gray')
axs[0, 1].set_title('Butterworth highpass filter, D_0 = 1%H, n = 4')
axs[1, 0].imshow(np_1_bw2i_4, cmap='gray')
axs[1, 0].set_title('Butterworth highpass filter, D_0 = 5%H, n = 4')
axs[1, 1].imshow(np_1_bw3i_4, cmap='gray')
axs[1, 1].set_title('Butterworth highpass filter, D_0 = 10%H, n = 4')
plt.show()
Наилучшие результаты опять достигаются при $D_0 \sim 5\%$.
Применим фильтр Гаусса для повышения резкости изображения.
np_1_gs1, np_1_gs1_fm = apply_gauss_HP(np_1_fft, .01)
np_1_gs2, np_1_gs2_fm = apply_gauss_HP(np_1_fft, .05)
np_1_gs3, np_1_gs3_fm = apply_gauss_HP(np_1_fft, .10)
np_1_gs1i = np.abs(get_ifft(np_1_gs1))
np_1_gs2i = np.abs(get_ifft(np_1_gs2))
np_1_gs3i = np.abs(get_ifft(np_1_gs3))
fig, axs = plt.subplots(2, 2, figsize=(15, 15), dpi=120)
axs[0, 0].imshow(np_1_sum, cmap='gray')
axs[0, 0].set_title('Original of the 1st image')
axs[0, 1].imshow(np_1_gs1i, cmap='gray')
axs[0, 1].set_title('Gaussian highpass filter, D_0 = 1%H')
axs[1, 0].imshow(np_1_gs2i, cmap='gray')
axs[1, 0].set_title('Gaussian highpass filter, D_0 = 5%H')
axs[1, 1].imshow(np_1_gs3i, cmap='gray')
axs[1, 1].set_title('Gaussian highpass filter, D_0 = 10%H')
plt.show()
Наилучшие результаты достигаются при $D_0 \sim 10\%$ -- например, становится хорошо читаемым слово susceptibility, также более четко видно букву $\omega$. Полученные изображения являются более гладкими, чем при использовании предыдущих двух фильтров. Гауссов фильтр дает хорошее качество фильтрации даже для маленьких объектов и тонких полос.